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We numerically study the superconductor-insulator phase transition in a model disordered 2D 
superconductor as a function of applied magnetic held. The calculation involves quantum Monte 
Carlo calculations of the (2 + l)D XY model in the presence of both disorder and magnetic Held. The 
XY coupling is assumed to have the form — J cos{9i — 9j — Aij), where Aij has a mean of zero and a 
standard deviation A Aij. In a real system, such a model would be approximately realized by a 2D 
array of small Josephson-coupled grains with slight spatial disorder and a uniform applied magnetic 
held. The different values AAij then corresponds to an applied held such that the average number of 
flux quanta per plaquette has various integer values A'^: larger TV corresponds to larger AAij. For any 
value of AAij, there appears to be a critical coupling constant Kc{AAij) — ^y[J/ {2U)]c, where U is 
the charging energy, above which the system is a Mott insulator; there is also a corresponding critical 
conductivity a*{AAij) at the transition. For AAij ~ oo, the order parameter of the transition is a 
renormalized coupling constant g. Using a numerical technique appropriate for disordered systems, 
we show that the transition at this value of AAij takes place from an insulating (I) phase to a Bose 
glass (BG) phase, and that the dynamical critical exponent characterizing this transition is 2; ~ 1.3. 
By contrast, z = 1 for this model at AAij = 0. We suggest that the superconductor to insulator 
transition is actually of this I to BG class at all nonzero AAij's, and we support this interpretation 
by both numerical evidence and an analytical argument based on the Harris criterion [A. B. Harris, 
J. Phys. C: Solid State Phys. 7, 1671 (1974)]. Kc is found to be a monotonically increasing function 
of AAij. For certain values of K, a disordered Josephson array may undergo a transition from an 
ordered, Bose glass phase to an insulator with increasing AAij. 



I. INTRODUCTION 



The superconductor-insulator (S-I) transition of thin 
two-dimensional (2D) superconducting films has been ex- 
tensively studied both theoretically Hi, [2, S II H, 1, 0, 
" "JTolJl^^ 13 El m, El and experimentally 
19l. 20r2ir22| for many years. The theoretical work 
can be broadly categorized into two groups: in one group, 
disorder is induced using a random chemical potential, 
while in the other, disorder is generated using a mag- 
netic field. Most previous work belongs to the former 



[J, i, i, i, 0, i, i, ITpJHJTiP whereas only a few 
belong to the latter [1 [llllllllf^ 

The present work is motivated primarily by several 
experiments in which an S-I transition is observed in a 
2D material as a function of applied transverse magnetic 
field. Such experiments have been reported in thin films 
of superconducting materials. They have also been car- 
ried out in some of the most anisotropic cuprate high-Tc 
superconductors; in such materials, individual copper ox- 
ide layers may conceivably behave like thin superconduct- 
ing films if they are well enough de coup led from the other 

layers [lllMlaEIMIlllMIMIMlMlMIMIMIMP- 

In both cases, the films seem to undergo a transition from 
S to I with increasing magnetic field. Furthermore, the 
transition appears to be controlled mainly by the film 
resistance R. Experiments suggest that, in contrast to 
some predictions, R does not have a universal value at 
the S-I transition 2^2M- view of these experiments, 
it seems useful to construct a simple model which con- 
tains disorder and also shows a field-driven transition. In 



the present paper we present such a model, and analyze 
its properties by a combination of numerical methods and 
scaling assumptions. 

Before describing our own approach, we briefly review 
some of the previous theoretical work in this area. An 
early numerical calculation was carried out by Cha et al. 
1;] at zero magnetic field (B = Q). These workers calcu- 
lated both analytically and numerically the zero-temper- 
ature (r — 0) universal conductivity a* of the 2D boson- 
Hubbard model without disorder at the S-I transition. 
They found, using numerical Monte Carlo (MC) simula- 
tions of a {2+l)B XY model, that a* = (0.285±0.02)crQ, 
where ctq = (2e)^//i is the quantum conductance. This 
result is close to the value obtained from an analytic 
large- iV expansion. They further studied this model un- 
der an applied transverse magnetic field using MC simu- 
lations, and found that a* was increased [2]. 

Fisher et al. 3 studied the T — phase diagrams 
and phase transitions of bosons with short-range repul- 
sive interactions moving in both periodic and random 
potentials. For the periodic case, they found the system 
exhibited two different phases, a superfluid and Mott in- 
sulator, and that the dynamic exponent z exactly equaled 
the spatial dimension d. They also derived certain zero 
temperature constraints on the correlation length expo- 
nent v and the order parameter correlation exponent 77, 
namely v > 2/d and 77 < 2 — d. In the presence of dis- 
order, they found that a "Bose glass" phase existed, and 
that the transition to a superfluid phase took place from 
the Bose glass phase, not directly from the Mott insula- 
tor. 
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Most previous studies in this area have been based on 
quantum Monte Carlo (QMC) simulations 0, i, i, S i, 
m E [H, H [ll ■ Some work has involved advanced 
QMC techniques, such as a QMC algorithm based on 
the exact duality transformation of the boson Hubbard 
model and a worm algorithm [14]. Other studies 
have used a stochastic series expansion method [ll|, [l3] 
and an exact diagonalization method p^ . Analytically, 
besides the large- iV expansion technique used in Ref. [l| , 
a coarse-graining approximation ^ has been adopted in 
some investigations. 

The numerical studies of the S-I transition have used 
a wide range of model Hamiltonians. Some workers have 
employed a 2D hard core boson model 0, |^, [ill, [H) 
while others used a 2D soft core boson Hamiltonian [l|, 
H, [Ij [33 • This model has been used to investigate the S-I 
transition at T = [1, [3, [1 [S [O, [H , as well as the su- 
perconductor-Bose glass (S-BG) phase transition [a, [S | , 
while some workers have investigated both [3, [l3j IT3LT14 j . 
In addition, Smakov and S0rensen [13] studied the S-I 
transition at finite temperature T using a similar model. 

A number of workers have also investigated more com- 
plex phase transitions, of which we mention just a few 
representative examples. Capriotti et al. [l^ studied a 
reentrant superconducting-to-normal (S-N) phase transi- 
tion using, as a model, a resistively shunted 2D Joseph- 
son junction array with normal Ohmic shunt resistors as 
the source of dissipation. Chakravarty et al. [11] also 
found a dissipation-induced phase transition in such an 
array, but did not study the possibility of reentrance. 
The reentrant S-N phase transition in Ref. [l^ was found 
to persist for moderate dissipation strength, but the su- 
perconducting phase was always found to be stabilized 
above a critical dissipation strength at sufficiently low T. 
Hebert et al. fll!] studied phase transitions between su- 
perfluid, checkerboard, and striped solid order, using two 
interactions — a nearest-neighbor (Vi) and next-nearest- 
neighbor (V2) repulsion — instead of a single parameter to 
describe the random chemical potential. They found that 
the model exhibited a superfluid to striped solid tran- 
sition at half filling; away from half filling, they found 
a first-order transition from superfluid to striped super- 
solid, as well as a continuous transition from striped su- 
persolid (superconducting) to striped solid (insulating). 
Schmid et al. 12] have studied a first order transition 
between a checkerboard solid and a superfluid phase at 
finite temperature. They also found that an unusual 
reentrant behavior in which ordering occurs with increas- 
ing temperature. As an effort to develop a more realis- 
tic model, several workers have included both short and 
long-range repulsive interactions between bosons [1, [13] , 
and some studies have included fluctuations in the am- 
plitude as well as the phase of the superconducting order 
parameter [7|. 

The T = S-I transition has been found to be charac- 
terized by universal behavior. Ref. [§] found, using QMC, 
that the dynamic exponent, the correlation length expo- 
nent, and the universal conductivity were z = 0.5 ±0.1, 



V = 2.2 ± 0.2, and Uc = (1.2 ± 0.2Wq, respectively. In 
the coarse-graining approximation [3] , the universal con- 
ductivity was found to be a* — (7r/8)(TQ, while the value 
a* = (0.45 ± 0.05)crQ was obtained at finite T using the 
stochastic series expansion with a geometric worm al- 
gorithm [l3]; in the latter work it was also found that 
(t/ctq scaled with uj /T at small frequencies uj and low T . 
With only short-range Coulomb interactions, the univer- 
sal conductivity at thephase transition was found to be 
a* = (0.14±0.03)crQ With long-range Coulomb 

interactions, this value increased to a* = (0.55 ± 0.1)crQ 
or cr* = (0.55 ± 0.06)crQ [l^. 

This critical behavior differs significantly from the 
T = S-BG transition. At this transition, the dynami- 
cal exponent and the universal conductivity were found 
to equal z = 1.95±0.25 and a* = (0.17±0.01)crQ, respec- 
tively [1]. Batrouni et al. [I found a* = (0.45 ± 0.07)crQ 
from QMC calculations and a* = (0.47 ± 0.08)(Tq from 
analysis of current-current correlation functions. 

At intermediate strength of disorder, Lee et al. [13] 
found that the dynamical and the correlation length crit- 
ical exponent were 1.35 ± 0.05 and v — 0.67 ± 0.03, re- 
spectively. They also found that a Mott insulator to su- 
perfluid transition occurred in the weak disorder regime 
while a Bose glass to superfluid transition took place in 
the strong disorder regime. More recently, Lee and Cha 
[1^ studied the quasiparticle energy gap near the quan- 
tum phase transition. They found that this gap vanished 
discontinuously at the transition for a weak disorder, im- 
plying a direct Mott insulator to superfluid transition, 
whereas this discontinuous jump disappeared for a strong 
disorder, supporting the intervention of Bose glass phase 
in this regime. 

Several workers have studied the S-I transition by ex- 
plicitly introducing a magnetic field, using various mod- 
els and experiments. For example, Nishiyama |15i] found 
that the 2D hard core boson model exhibited a field- 
tuned localization transition at a certain critical mag- 
netic field and that the critical DC conductivity was 
substantially larger than that at zero magnetic field. In 
his work, the critical conductivity was found to be non- 
universal but instead increased with increasing magnetic 
field. Besides the experiments mentioned earlier, Sam- 
bandamurthy et al. [l9[ found, from studies of thin amor- 
phous InO films near the S-I transition, that the resistiv- 
ity followed a power law dependence on the magnetic field 
in both the superconducting and the insulating phases. 

In most of the above calculations, the QMC approach 
is based on a mapping between a d-dimensional quantum 
mechanical system and a (d-l-l)-dimensional classical sys - 
tem with the imaginary time as an extra dimension [391 ]. 
This mapping works because calculating the thermody- 
namic variables of the quantum system is equivalent to 
calculating the transition amplitudes of the classical sys- 
tem when they evolve in the imaginary time. The imagi- 
nary time interval is fixed by the temperature of the sys- 
tem. The net transition amplitude between two states of 
the system can then be obtained by a summation over the 
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amplitudes of all possible paths between them according 
to the prescription of Feynman These paths are 

the states of the system at each intermediate time step. 
Therefore, the path-integral description of the quantum 
system can be interpreted using the statistical mechan- 
ics of the {d + l)-dimensional classical system held at a 
fictitious temperature which measures zero-point fluctu- 
ations in the quantum system. 

In order for a boson system to have a superconduc- 
tor-insulator transition, the bosons must have an on-site 
repulsive interaction, i.e., a "charging energy." Other- 
wise the boson system would usually undergo Bose-Ein- 
stein condensation at zero temperature. The charging 
energy induces zero-point fluctuations of the phases and 
disorders the system. On the other hand the Josephson 
or XY coupling favors coherent ordering of the phases, 
which causes the onset of superconductivity. Therefore, 
the competition between the charging energy and XY 
coupling is responsible for the superconductor-insulator 
transition, which typically occurs at a critical value of the 
ratio of the strengths of these two energies. In addition, 
if a disorder is added to the system, the system may also 
undergo a transition to a phase other than a Mott insu- 
lator, depending on the strength of the disorder. This 
additional phase is known to be a Bose glass phase. 

In this work, we study the zero-temperature quan- 
tum phase transitions of 2D model superconducting films 
in an applied magnetic field. Our model includes both 
charging energy and Josephson coupling, and thus allows 
for an S-I transition. In our approach, the applied mag- 
netic field is described by a root-mean-square (rms) fluc- 
tuation AAij which describes the randomness in the flux 
per plaquette. This randomness leads to the occurrence 
of a Bose glass phase at large AAij. As explained fur- 
ther below, our model corresponds well to a 2D Josephson 
junction array with weak disorder in the plaquette areas, 
studied at an applied uniform magnetic field correspond- 
ing to integer number Ny, on average, of flux quanta per 
plaquette. The quantity AAij is proportional to the rms 
disorder in the flux per plaquette, and is proportional to 
Ny. Thus, our model gives rise to an S-I transition in the 
array with increasing Ny (or increasing magnetic field) . 

The remainder of this paper is organized as follows. 
Sec. In] presents the formalism. In this section, we give 
the model boson Hubbard model, and describe its con- 
version to a (2 -I- 1)D XY model, which we treat using 
path-integral Monte Carlo calculations. We also describe 
the finite-size scaling methods for obtaining the critical 
coupling constants and universal conductivities at the 
transition. Finally, this section describes the nature of 
the renormalized coupling constant used to study the be- 
havior of the system in the fully random case. Sec. IIIII 
presents our numerical results, using these approaches. 
We discuss our results and present our conclusions in 
SecHVl 



II. FORMALISM 



Model Hamiltonian 



Our goal is to examine the superconductor-insulator 
transition in a disordered 2D system in a magnetic field 
at very low temperature T. Thus, a useful model for this 
transition would include three features: (i) a competition 
between a Coulomb energy and an energy describing the 
hopping of Cooper pairs; (ii) disorder; and (iii) a mag- 
netic field. In particular, we hope that this model will 
exhibit, for suitable parameters, a transition from super- 
conductor to insulator as the magnetic field is increased. 
While there are a wide range of models which could in- 
corporate these features, we choose to consider a model 
Hamiltonian appropriate to a 2D Josephson junction ar- 
ray: 



J^cos(0i - 6 J - Ai 



(1) 



Here Uj is the operator representing the number of 
Cooper pairs on a site j, J is the Josephson energy cou- 
pling sites i and j , 9i is the phase of the order parameter 
on the ith site, Aij — (27r/<I>o) A • dl is a magnetic 
phase factor, $0 — hc/2e is the flux quantum, and A 
is the vector potential. In this picture, each site can be 
thought of as a superconducting grain. 

For calculational convenience, we choose to take the 
sites j to lie on a regular 2D lattice (a square lattice in 
our calculations), with Josephson coupling only between 
nearest neighbors. Thus, the disorder in our model is 
incorporated via the magnetic phase factors Aij, as ex- 
plained further below. Our Hamiltonian is identical to 
that of Cha et al. [ij except that we consider the special 
case that the chemical potential fii for Cooper pairs on 
the ith grain is an integer, and we choose the ^^'s to be 
random. 

The first term in Eq. ^ is the charging energy. We 
consider only a diagonal charging energy and also assume 
all grains to be of the same size, so that U is independent 
of j . Since the charging energy Ecj of a grain carrying 
charge Qj with capacitance C is Ecj = Q|/(2C), 



U 



(2e)2 _ 2e^ 
2C " IT 



(2) 



We also know that Qj — CVj , where Vj is the voltage of 
grain j relative to ground; so 



2(2e) 



2 3 ' 



(3) 



where we have used the Josephson relation, Vj = 
{fi/2e)0j. Finally, we can express C in terms of U us- 
ing Eq. ([2]), with the result 



(4) 
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Combining all these relations, we obtain 
AU 



(5) 



Since we have taken the grains to lie on a lattice, we 
need to choose the A^-'s in a way which incorporates 
randomness. Thus, we make the simplifying assumption 
that the phase factor Aij of each bond in the plane is an 
independent Gaussian random variable with a mean of 
zero and a standard deviation AA,-: 



1 



exp 



^1 



2(A^,,)2 



Since the sum of the phase factors around the four bonds 
of a plaquette is 27r/$o times the flux through that pla- 
quette, this choice will cause the flux through the pla- 
quette also to be a random variable. However, the fluxes 
through nearest-neighbor plaquettes will be correlated. 

Although this model may seem artificial, it should 
closely resemble a real, physically achievable system. 
Specifically, consider a spatially random distribution of 
grains in 2D in a uniform magnetic field. Suppose that 
the grain positions deviate slightly (but randomly) from 
the sites of a square lattice. Then the areas of the square 
plaquettes have a random distribution, and consequently, 
the flux through each plaquette also varies randomly 
about its mean. If the average flux $ per plaquette is 
^ = f^o, then the root-mean-square deviation of the 
flux, Acf> cx / as in Fig. [H 




exactly like the array at / = 0, because the Hamiltonian 
would then be perfectly periodic in / [ill, 113, Ell ■ With 
nonzero disorder, only the rms deviation from integer /, 
i.e.. A/ = A$/$o, is physically relevant. This deviation 
increases linearly with /. 

In short, our model Hamiltonian is approximately real- 
ized by a 2D Josephson junction array on a square lattice, 
in which the grains deviate randomly in position from 
their lattice sites, placed in a transverse magnetic field 
with an average flux per plaquette /$0i with integer /. 
A larger A Aij corresponds to a larger /. The models are 
not equivalent, even if the random position deviations are 
specified by Gaussian variables, because we assume the 
(6) AAij 's for different bonds are uncorrelated, whereas they 
would be correlated in the positionally disordered case. 
However, this difference should have little effect in prac- 
tice and we have confirmed this in the limit of large / (see 
below). For non- integer /, the disordered 2D Josephson 
array has well-known oscillatory properties as a function 
of / which are not described by our model as formu- 
lated above. This disordered Josephson array is not an 
entirely realistic model of a superconducting film which 
undergoes a field-driven superconductor-insulator transi- 
tion because our model involves an underlying lattice of 
grains. Nonetheless, we may hope that some of the prop- 
erties of our model resemble those seen in experimentally 
studied materials. 

A related method of including random flux has previ- 
ously been used by Huse and Seung [44] as a model for 
a three-dimensional (3D) "gauge glass." These workers 
considered only A Aij ~ oo, and studied a 3D classical 
model (U = 0) rather than the 2D quantum case consid- 
ered here. 



B. Path Integral Formulation 

We can now use the model Hamiltonian (O to obtain 
the action in the form of a standard integral over imagi- 
nary time. The action S may be written 



FIG. 1: Sketch of a 2 x 2 group of plaquettes in a square lattice 
of Josephson-coupled grains, in which each grain is randomly 
displaced by a small amount from its nominal lattice site. In 
a uniform transverse magnetic field, if the average flux "1> per 
plaquette is <1> = f^o, then the root-mean-square deviation 
of the flux from its mean value is also proportional to /. 

Now consider the special case of integer /. In the ab- 
sence of disorder, the 2D array at integer / should behave 



S 

h 



1 



Cdr, 



(7) 



where C is the Lagrangian, given by 



C 



4C7 ^ 



dr 



- J^cos[0,(t)-0,(t)-A,,(t)]. 

(8) 

The partition function is now given by a path in- 
tegral of exp{~-S/'h) over all possible paths described 
by the variables Oi(T) in imaginary time r, integrated 
from T = to T ^ l3h, where f3 = l/{kBT). This 
path integral can be reduced to the partition function of 
an anisotropic classical XY model in three dimensions. 
Here, by "anisotropic" we mean that the coupling con- 
stants K and Kr in the xy plane and r direction are 
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different. To make the mapping, we first write 



1. Helicity Modulus 



dr 



A0, 
At 



2 - 2 cos A9, 
(Ar)2 



(9) 



where At is the width of the time shce, A9i = 0i{T + 
At) — 9i (r) , and we have used the expansion of cos Ad to 
second order in the smaU quantity AO. This expansion 
is accurate when At is sufficiently small. 

Neglecting the constant term in this expansion, we fi- 
nally obtain 

I = -if,^COs[0,(T)-0,(T + AT)] 

-KJ2co8[9.,{t) - e,{T) - A,,(t)]. (10) 

Here the sums run over all bonds in the t direction and 
in the xy plane, respectively. 

In order to obtain the values of the coupling constants 
K and K^^ we assume that we have broken up the time 
integral into M time slices, each of width j3Ti/M . Then 
the coupling constant in the xy direction is just 



K 



(11) 



The coupling constant in the t direction is given by 
2 I h M 



1 



4Uh^ (At)2 2U At 2/3C/' 



(12) 



we have used At = (3fi/M and included the extra factor 
of 2 in the numerator because cos A0 ~ l — (A0)^/2. 

Given K and Kr, the partition function is obtained 
from this anisotropic 3D XY classical Hamiltonian with 
coupling constants K and Kt. Any desired equilibrium 
quantity can, in principle, be computed by averaging over 
all configurations using standard classical Monte Carlo 
techniques. Within any given realization of the disorder, 
the Aij 's are chosen at random from the Gaussian distri- 
bution within the xy plane, as described above, but the 
Aij 's for a given bond in the xy plane are independent of 
T, i.e., they are the same for all time slices. In principle, 
for any given /3, this should be done taking the limit as 
M — > oo. In practice, of course, the size of the sample is 
limited by considerations of computer time. 



C. Evaluation of Specific Properties Using Path 
Integral Formulation 

The time-slice formulation of the partition function al- 
lows various properties to be evaluated using standard 
classical Monte Carlo techniques. We now review how 
this may be done for the helicity modulus (or superfluid 
density) and the electrical conductivity. Similar formu- 
lations have been given in Refs. [11, |3, IS IH, [45|, |46| for 
different but related models. 



For a frustrated classical XY system in d dimensions, 
the helicity modulus tensor is a d x d matrix which 
is a measure of the phase stiffness. It is defined as the 
second derivative of the free energy with respect to an 
infinitesimal phase twist, and may be written 



7a/3 



1 



TV dA'dA 



(13) 



where N is the number of sites in the system and A' is a 
fictitious vector potential added to the Hamiltonian (in 
addition to the vector potential A already included in 
the Hamiltonian [1^). In explicit form, this derivative 
takes the following form for the diagonal elements (see, 
e.g., Ref. (13): 



aa — AT \ ^ ^ '^^j ^^^i^i ^ii)(^ii ' ^a) 




,-d,-A,,){e,j-ec.)) .(14) 



Here is a unit vector from the zth to the jth site, and 
Cq is a unit vector in the a direction. The triangular 
brackets denote an average in the canonical ensemble. 

If this expression is applied to the time-slice represen- 
tation of the quantum-mechanical Hamiltonian, the cou- 
pling constants Jy will be different in the xy plane and 
in the t direction. For the time-slice calculation, we have 
to be careful in order to obtain a result which is well-be- 
haved in the limit M 00, where M is the number of 
time slices. The correct expression in this case is 



9, -A,) 



9 J Aij ) 



A/}) . (15) 



Here we are assuming that there are N^Ny superconduct- 
ing grains in our 2D lattice and M time slices. The sums 
run over all distinct bonds in the x direction; there are 
NxNyM of these bonds {N^Ny per time slice). A similar 
expression holds for "fyy. The in-plane coupling constant 
is taken to be J/M because there are M time slices. 

From the above form, we can see why the expression 
behaves correctly in the limit M ^ 00. Each of the 
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two sums contains N^NyM terms in it, but the ensemble 
average consists of Af identical terms, one for each layer. 
Therefore, the first sum in Eq. (fTSjl . for example, should 
take the form 




(16) 



where the sum on the right hand side runs only over the 
phases in a single layer. The right-hand side is evidently 
independent of M in the limit M oo. A similar ar- 
gument can be used to show that the second part of the 
expression for 'y^x also approaches a well-behaved 
limit as M —^ oo. Our numerical results confirm this 
behavior. 

As an illustration, we write down an expression for 
in the limit T ^ in the unfrustrated case (A Ay — 0). 
First, we multiply expression (fT5)) by f3/M to obtain 



Now it is known from previous Monte Carlo studies 
m, m, [13, [m, m that the unfmstrated 3D XY model 
on a simple cubic lattice has an ordered phase if i^T > 
Kc - 1/2.21 - 0.452. Therefore, ^{K) vanishes if K < 
Kc and is positive for K > Kc- Translating this result 
to the 2D quantum XY model on a square lattice, we 
see that there is a superconductor-insulator transition at 
J/{2U) = (0.452)2 = 0.204. 

For reference we give the connection between our for- 
mulation of the helicity modulus and the calculation of 
Cha et al. [l| Rather than the helicity modulus, these 
workers calculate the quantity /9(0), which is related to 
the superfluid density ps by 



P(0) 



k^T 



(20) 



and to the components of the helicity modulus tensor by 
p{0) = = K^-y^x + 7yy)/2, where ^xx = 7yy for the 
present model, which is isotropic in the xy plane. In our 
notation, p{0) is given by 



Plx 



M 




N^NyM 



N^NyM 



N^NyM 



^ Xsin( 

\(ij)\\x 



0, - A; 



Aij) 



(17) 



where K — fiJ/M . The corresponding coupling constant 
in the t direction is Kr = M/ {2(311). 

Since we are interested in the limit /? » 1, we choose 
M so that K = Kr. This condition is equivalent to 



M 



1 



V2JU 



(18) 



Hence, we get 



y/2JU N^NyM 
1 



^ IkY, cos(a,-0,-A,,) 



N^NyM 



(u>ll* 
(«i>ll* 



N^NyM 



(19) 



where K = /3J/A/ = ^J/(2C/). Since K ^ Kr, the 
right-hand side of this equation represents a dimension- 
less helicity modulus 7 for a classical unfrustrated 
isotropic 3D XY model on a simple cubic lattice, which 
is a function of a single dimensionless coupling constant 
K. 



P{0) 



JK 



2NxNyM 



Y cos{0i - 0j - Aij) 

{ij)\\x I 



JK' 



2N.,NyM 



Y sin(6ii - 0.j - A^ 

(«i>ll* 



Y sin(6l, - Oj - A, 



+ 



JK' 



2NxNyM 



{ij)\\x I 

Y smiOi-Oj- Aij)\ 
\{i3)\\y I 



2. Specific Heat 



(21) 



For the specific heat Cy, we used the fluctuation-dis- 
sipation theorem given by 



Cv = 



NkeT' 



(22) 



where iV is the total number of sites in the lattice, Ti' is 
the Hamiltonian in Eq. and (• • ■ ) denotes an ensem- 
ble average. 
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3. Conductivity 

The conductivity in the low-frequency hmit can also 
be obtained from the time-slice Monte Carlo approach 
as [H 

pik) 



cr(0) = 27r(7Q lim 



k^O k 



(23) 



where p(k) is proportional to the superfluid density at 
frequency k and is given by 



Pik) = 



JK 



2N^NyM 




cos{9i - 9j - A, 

\{ij)\\x 



A-ij) 



sm{0i - Oj - Aij)e~ 



2N^N„M 



J2 sh 

(ii>IIS,x 



9j - A,j)e 



X ^ sin(0i — — Ai^)e 

(u)llfi.x 



2Na,NyM 



sin(a,-(?,-A,,)e-*-^ 




A,j)e 



(24) 



In the limit of very small fc, we expect that p{k) will 
remain finite in the superconducting phase and vanish in 
the insulating phase. Thus, cr(0) will become infinite in 
the superconducting phase but vanish in the insulating 
phase. Precisely at the critical value Kc, cr(0) will become 
finite with a universal value, as already obtained by other 
workers for related models. 



III. QUANTUM MONTE CARLO RESULTS 

A. Numerical Procedure 

In our quantum Monte Carlo calculations, we use the 
standard Metropolis algorithm with periodic boundary 



conditions in both the spatial directions and the imag- 
inary time direction. We usually start with a random 
configuration of phases at K — 0.4, then increase K up 
to K ~ 0.7 in steps of 0.005. This procedure corresponds 
to lowering the temperature T since K (x 1/T. At each 
K, we take 40000 MC steps per site through the entire 
lattice to equilibrate the system, after which we take an 
additional 50000 MC steps to calculate the thermody- 
namic variables of interest. For a lattice size of 6'^, we 
use ten times as many MC steps as these for both equili- 
bration and averaging, and for a lattice size of 8^, we use 
twice as many steps. 

For the phases 9i of the order parameter on each site, 
we use the 360-state clock model instead of a continuous 
angle between and 27r since it allows us to cover the 
entire range of angles with fewer trials. Therefore, the 
allowable angles are 0°, 1°, 2°,..., 359°. It has been shown 
numerically that these discrete phase angles give results 
indistinguishable from the continuous ones provided that 
n > 20 [53|. However, we select Aij from a continuous 
distribution in all our calculations. 

For the partially random and completely random Aij , 
we averaged over 100 different realizations of AAij to 
calculate the helicity modulus 7 and the specific heat Cy . 
These calculations were so time-consuming that we could 
go just up to 20 X 20 X 20 lattice size. For this reason, 
we chose to carry out simulations only over four different 
AAij for the partially random case: AAij — 1/2, l/\/2, 
(l-|-l/-\/2)/2 « 0.854, and 1. Each realization is specified 
by a different random number seed. 



B. Finite-Size Scaling for 7(0) 

In general, if there is a continuous phase transition 
as a function of some parameter, such as the coupling 
constant if, the critical behavior near the transition can 
be analyzed by carrying out a finite-size scaling analysis 
of various calculated quantities. For example, the zero- 
frequency helicity modulus 7(0) is expected to satisfy 



7(0) = 



1 



(25) 



where d is the spatial dimensionality, z is the dynamic 
exponent, 7 is a scaling function, v is the critical ex- 
ponent for the correlation length ^, (5 = {K — Kc)/Kc: 
Kc is the critical value of the coupling constant, and Lr 
is the thickness in the imaginary time direction. For 
our present system, d = 2, so the right-hand side is 
L-^^{L^/''6,Lr/L^). If we define ^{L^/''S,Lr/L^) = 
{L^ lLr)G{L'^/"5,Lr/L^), then this scaling relation be- 
comes 



L.7(0) = g('lV''<5,^ 



(26) 



If our computational box has N^j. x iV^ x M sites, with 
Nx — Ny, this scaling relation may be equivalently writ- 
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ten as 

KMj{Q) = G ^) . (27) 

As has been noted by other workers 
[13, HEHH, Kc can now be found, if 7 is a suitable order 
parameter, by pfotting KMj{0) as a function of S for var- 
ious cell sizes, all with aspect ratios satisfying M = ciVJ, 
and finding the point where these all cross, which cor- 
responds to 5 = 0. Unfortunately, this method requires 
knowing the value of z in advance. For many such quan- 
tum phase transitions, z may not be known. Thus, one 
should carry out this calculation for all plausible values 
of z and find out which value leads to a satisfactory cross- 
ing. This procedure is prohibitively demanding numer- 
ically. We have, therefore, initially attempted to carry 
out scaling using z — 1, the value which is known to 
be correct at AAij =0. If this value were correct also 
at AAij 0, it would suggest that the superconduc- 
ting-insulating transition at finite AAij is in the same 
universality class as the zero-field transition. In prac- 
tice, we find that z = I never gives perfect scaling at 
nonzero AAij and the scaling fit becomes progressively 
worse as AAij increases. For large AAij in particular, 
the fit clearly fails, and we find for these values that 7(0) 
never converges to a nonzero value. At such large AAij, 
we suggest that this regime corresponds to a Bose glass 
(as further discussed below), and carry out a different 
kind of scaling calculation to obtain the actual value of 
z of this phase at AAij = 00. We also give arguments 
suggesting that, in fact, this Bose glass phase is actually 
the ordered phase for all nonzero values of AAij ■ 

Operationally, we implement the hypothesis that z — 1 
by taking ^ Ny ^ M. With this choice, Eq. ((27l) 
becomes 

KMj{0) = GiM^'^S, 1) (28) 

when d — 2. 



C. Zero Magnetic Field 

As a check of our method, we have calculated Kc for 
the case of zero magnetic field {Aij = 0), using the above 
numerical approach. When there is no magnetic field, 
we get Kc = 0.4543 ± 0.0011 using a finite-size scaling 
analysis of the helicity modulus 7 as shown in Fig.[2l This 
value is very close to Kr = 0.4539 ± 0.0013 by the series 
expansion as in Ref. [53, which is also used in Ref. [H. 
This result confirms the validity of our numerical codes. 
Using our value of we can also obtain the universal 
conductivity cr* jaq = 0.282 ± 0.005. This result is also 
very close to the value a* jaq — 0.285 ± 0.02, obtained 
in Ref. f?]. 



■ 6x6x6' ' 
■r 8x8x8 

- * 1 ux 1 Ox U) 

♦ 12x12x12 

• 16x16x16 

- • 20x20x20 











0.43 0.44 0.45 0.46 0.47 0.48 



K 



FIG. 2: Plot of 7l/T as a function of K for various A^^, x 
Ny X M lattices for Aij = 0. In this and all subsequent 
figures, unless otherwise specified, we use N:^ ~ Ny = M. 
The phase transition occurs where the curves of different M 
cross. The crossing point yields Kc = 0.4543 ± 0.0011. The 
use of Nx = Ny — M is equivalent to assuming that the 
dynamic exponent z = 1, as discussed in the text. 

D. Finite AAij 

Figs. El^a), [3I^b), and m show the helicity modulus 7, 
the specific heat Cy, and the finite-size scaling behavior 
of 7 as a function of coupling constant K for several 
lattice sizes when AAij = 1/2. When K > 0.55, 7 and 
Cv appear to be nearly lattice size i ndep endent. The 
error bars from the jackknife method [58| are shown in 
Fig. EKa), but they are smaller than the symbol sizes. 
The lines are cubic spline fits to the data in Fig. ^h). 
The apparent crossing point in Fig. [?] yields the critical 
coupling constant Kc — 0.491 ±0.001, which is very close 
to the peak of Cy in Fig. [3Kb). 




FIG. 3: (a) Helicity modulus 7 and (b) specific heat Cv, 
plotted as functions of coupling constant K for several lattice 
sizes when AAij ~ 1/2. The error bars in (a), as obtained 
from the jackknife method, are smaller than the symbol sizes. 
The lines in (b) are cubic spline fits to the Monte Carlo data. 

The corresponding results for AAij = l/\/2 are shown 
in Figs. EKa), E^b), and [51 Compared to Fig. [SJa), 7 
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FIG. 4: Finite-size scaling behavior of the data in Fig. ISja), 
using z = 1. The apparent crossing point yields Kc = 0.491 ± 
0.001. 



FIG. 6: Same as Fig. |4l except that the data are from Fig. 
\^a.). The apparent crossing point yields Kc = 0.533 ± 0.001, 
as indicated by the vertical dashed line. 



shows more lattice size dependence when K > Kc in 
Fig. EJa). The apparent crossing point in Fig. [H] yields 
Kc — 0.533 ± 0.001. This Kc is also very close to the 
peak in Cy as in Fig.ISjb). 
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(b) 





0.40 0.45 0.50 0.55 0.60 0.65 0.70 

K 

FIG. 5: Same as Fig. [3 except that Aylij = 1/^/2. 

The results for 7, Cy, and when AA^ — 0.854 
are shown in Figs. El^a), El^b), and [H In this case, the 
lattice size-dependence of 7 when K > Kc in Fig. [7]^ a) 
becomes more conspicuous than that of Fig. [5l^a). The 
apparent crossing point for the different sizes in Fig. [8] 
is less clearly defined than in the previous examples, but 
yields Kc = 0.585 ±0.004. This Kc is slightly larger than 
the value of K at the maximum of the broad peak in 
Cy, as in Fig. [Tfb). There are some fluctuations of 7L7- 
around K — 0.70 for the lattice size of 12'^ and larger 
fiuctuations above Kc for the lattice sizes of 16'^ and 20^^ 
in Fig. m 

The fact that the apparent crossing point in Fig. [8] 
is even less clear than those for smaller values of A^l^ 
suggests that z ^ 1. We present a more likely scenario 
for this and other values of AAy in the discussion section 
below. 

As a final calculation for partially random Aij , we use 
A^ij = 1.0. The corresponding three thermodynamic 




0.40 0.45 0.50 0.55 0.60 0.65 0.70 



K 



FIG. 7: Same as Fig. |3] except that AAij = 0.854. 




0.40 0.45 0.50 0.55 0.60 0.65 0.70 



FIG. 8: Same as Fig. 14] except that the data are from Fig. 
[3a). The apparent crossing point yields Kc ~ 0.585 ± 0.004. 



variables 7, Cy, and jLr are shown in Figs. [SI a), M)^), 
and [TOl respectively. The lattice size-dependence of 7 in 
Fig-IHa) becomes far more conspicuous than the previous 
two examples. The peak of Cy is very broad, as shown 
in Fig. [Hb). In addition, there is nothing like a clear 
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crossing point of ^Lr for different sizes in Fig.[TOl We 
interpret this result to mean that the hehcity modulus 7 
does not play the role of an order parameter and that the 
transition is not a superconductor-to- insulator transition 
of the same character as at /S.Aij — 0. Furthermore, 
there are strong fluctuations of ^Lt as a function of Nx 
when K > 0.64 for most lattice sizes, as can be seen in 
Fig. [ini We believe that, for this value (and, in fact, at 
all nonzero values) of AAy, this is a transition from a 
Bose glass to a Mott insulator. 



other work fE^ , in the context of a different model. As 
in Figs.inja) andfTUl 7 is strongly lattice-size dependent 
and there exists no value of K at which the curves of 
Lt^{K) for different N.^ all cross (we do not show a plot 
exhibiting this lack of crossing). All these results indicate 
that we need a different order parameter to describe the 
phase transition. Besides these results, we find that the 
peak in Cv{K) is even broader than that in Fig. HKb). 
Moreover, the peak of Cv shifts towards a larger value 
of K as AAi, increases. 



0.12 
0.10 

0.08 
g 0.06 
="0.04 
0.02 
0.00 
1.6 
1.4 
•^1-2 
1.0 



6x6x6 

8x8x8 
• 10x10x10 
-• 12x12x12 
— 16x16x16 

20x20x20 



(a) 



■ ;.••:;••!':; 




0.03 
0.02 
0.01 
0.00 

-0.01 
0.04 
0.03 
0.02 
0.01 
0.00 

-0.01 



* 20X20X20 _ •.*!.#•.;, 

ii*H»<iHfHiiMff»HiHil{tllt|ft||f|{|"*:'** • i, , 



6x6x6 
8x8x8 
• 10x10x10 
12x12x12 
16x16x16 
20x20x20 



(a) 



6x6x6 

^" 8x8x8 

• 10x10x10 
12x12x12 
16x16x16 
20x20x20 



(b) 



0.45 0.50 0.55 0.60 0.65 O.70 
K 



FIG. 9: Same as Fig. [3] except that AAij = 1.0. 
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FIG. 10: Same as Fig. [l] except that the data are from Fig. 
|9l[a). In this case, the plots of Lt'^{K) for different A^^, do not 
cross, suggesting that the helicity modulus 7 is no longer a 
suitable order parameter at AAij = 1.0. 

Finally, we have considered the case of a fully random 
Aij^ A Ay — 00. We implemented this by choosing AAij 
randomly between and 27r. The helicity modulus 7 and 
the specific heat Cy for this case are shown in Figs.fTlTa') 
and [TWa) , respectively. The magnitude of 7 becomes 
much smaller than those of previous cases, so that the 
error bars are easily visible on the scale of the plot. The 
helicity modulus even seems to have negative values for 
certain values of K, depending on the lattice size. Such 
negative values and fluctuations of the helicity modulus 
in a disordered superconductor were already reported in 



FIG. 11: (a) The helicity modulus 7 as a function of coupling 
constant K for several lattice sizes and AAij = 00. (b) Same 
as (a) except that we assume a uniform transverse magnetic 
field with frustration f — 20 and disorder in the grain po- 
sitions with a uniformly distributed random displacement of 
each site, as described in the text. In both (a) and (b) we find 
a negative 7 for certain values of K, depending on the lattice 
size. The similarity of (b) and (a) is evidence that these two 
models give very similar results. 
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FIG. 12: Same as Fig. [TT] but for the specific heat Cv- The 
lines are cubic spline fits to the data. 

As a comparison to the fully random AAij case, we 
have also considered a model of "positionally disordered 
sites" in a strong uniform transverse magnetic field B = 
Bz, similar to a model considered in Ref. The po- 

sition coordinates {xi,yi) of each site are assumed uni- 
formly and independently distributed between —A and 
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A with respect to the position 
have in the ordered lattice, i.e., 



'xio,yio) the site would 



\xi -Xio\ < A, 
lUi - y«| < A. 



(29) 



In our calculations, we have chosen A = a/4, where a is 
the lattice constant of the unperturbed lattice. Thus Aij 
has the form 



_2Tr Xi + 



-{Vj - Vi) 



(30) 



for nearest-neighbor sites i and j. In order to consider 
a strong field, we choose / = Ba^/^o = 20. The results 
for this system of positionally disordered sites are shown 
in Figs. fllTb) and [TWb). They are qualitatively similar 
to those with AA^j — oo, and even quantitatively simi- 
lar for Cv- We conclude that the model of positionally 
disordered sites is nearly equivalent to that with random 
Aij, at least for A Aij — oo. 

At large values of A Aij, our results suggest that the 
transition occurs between a Mott insulator and a Bose 
glass rather than a conventional superconductor. Since 
a new order parameter is demanded to study this tran- 
sition, we use the "renormalized coupling constant" g as 



Ref. 



Using the same Aij for each realization, two 



replicas of phase 6j are simulated with different initial 
conditions and updated using different random numbers. 
Their overlap is calculated from the quantity 



Vexp[z(0W-0f )], 



(31) 



where and Oy are the phases at site j in the two 
replicas. Given q, the renormalized coupling constant g 
is defined as 



n(2) 



5 = 2 



191- 



(32) 



where (• • • ) denotes the thermal average while [• • •] de- 
notes an average over many realizations of Aij . Figure [T3l 
shows this 5 as a function of coupling constant K for 
several lattice sizes when A Aij — oo. From the crossing 
point for different sizes, we obtain Kc = 0.630 ± 0.005. 
Unlike the results in Ref. [31, g still has a size depen- 
dence when K > Kc- 

Figs. [TlHTSl strongly suggest that AAij — oo corre- 
sponds to a Bose glass transition, rather than a con- 
ventional superconducting transition. Hence, we expect 
z ^ 1. In order to allow for z ^ 1, we have carried 
out additional calculations, using a method suggested by 
Guo et al. [g^l and by Rieger and Young t6ll|. Follow- 
ing the procedure of these authors, we first calculate g 
as a function of the time dimension L^- for various sizes 
L and several temperatures T. Since the proper scaling 
behavior of g is not expected to depend on the anisotropy 
of the coupling constants, we assume the same coupling 
constant J = 1 in both the space and imaginary time 




FIG. 13: The renormalized coupling constant g [Eq. (|32|l ] as a 
function of coupling constant K for several lattice sizes when 



oo. The crossing point yields = 0.630 ± 0.005. 



The lines are cubic spline fits to the data. 



directions. For each T and L, g has a maximum value 
as a function of Lr- According to Refs. [6^ and [6l|, the 
true Tc is the temperature such that this maximum value, 
5max, is independent of L. Once Tc is determined by this 
procedure, the correct z is that value which causes a plot 
of g{Tc, Lr/L^) versus the scaling variable Lr/L^ to be 
independent of L. 

Following this prescription, we have calculated 
g{Lr, L, T) as a function of Lt for various values of T and 
a lattice of size L x L x Lr, assuming that J ^ Jr = I- 
We find that ^max is most nearly independent of L when 
T — l.GlJ/fcs. To illustrate this independence, we plot 
g{Lr,T = 1.61J/A;b) for several choices of L in Fig. [TH 
At this temperature, for all L studied, g{Lr,L,T) has a 
maximum of around 0.38 when plotted against Lr- 
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FIG. 14: Plot of g[Lr, L, T) versus Lr for several values of L, 
as given in the legend, for T = 1.61J/fcfl and AAij = oo. In 
these calculations, the coupling constants J and Jr are each 
taken to be unity. For each L, each calculation represents an 
average over 100 realizations of the disorder. The maximum 
values gina.x{L) are nearly independent of L. 

Given T^, we obtain z by plotting g(Ti., Lr/L^) as a 
function of the scaling variable Lr / L' for various values 



12 



of L. The correct value of z is the one which causes these 
curves to be most nearly independent of L. We have 
made such plots for various values of z aiTc — 1.61J//cb, 
and find that this collapse of the numerical data is most 
nearly obtained for z — 1.3, with an uncertainty of about 
±0.1. The resulting scaling fit is shown in Fig.[T5]for z — 
1.3. The fit is very good, suggesting that (i) the transition 
for AAij = oo is indeed a Bose glass transition, and (ii) 
the critical exponent z at the transition is z ^ 1.3 ± 0.1. 
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FIG. 15; Plot of g{Lr,L,T) versus logjg(L.r/i^) for several 
values of L, as given in the legend, for T — 1.61J/fc_B, z = 1.3, 
and AAij = oo. In these calculations, the coupling constants 
J and Jt are each taken to be unity. For each L, each point 
represents an average over 100 realizations of the disorder. For 
this choice of z, the results for different values of L collapse 
very well onto a single plot. The corresponding plots for z — 

1.2 and z = 1.4 produce only slightly inferior collapses. We 
conclude that the correct value of z for this transition is z ~ 

1.3 ±0.1. 
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FIG. 16: Same as Fig. [Til except for AA^j = 1.0 and T = 
1.70 J/fcs. 
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FIG. 17: Same as Fig. [15] except for AA^j = 1.0 and T = 
1. 70 J/fcs. 



We have carried out a similar series of calculations at 
AAij = 1-0. For this choice, the best glass scaling fits 
are reasonable, but not so good as for AAij = oo. They 
are shown in Figs. [TBI and [T71 for T = 1.70 J/ ks, which 
is our best estimate for the glass transition temperature 
of this model at AAij = 1-0. Our conclusion is that, 
for AAij = 1-0, the sizes we can achieve [L ^ Lr ^ 12) 
are simply not large enough to reveal the excellent scaling 
behavior which is expected for a sufficiently large sample. 
We discuss below a possible explanation why AAij = 1-0 
requires a larger sample size than AAij — oo. 

With all the KcS we have collected so far, we can 
plot \/Kc as a function of AAij. This is shown in Fig. 
rrSl Since Kc — ^J[J / {2U)]c and since 1/Kc decreases 
as AAij increases, these results mean that [J/{2U)]c in- 
creases with increasing AAij. Therefore, there exist cer- 
tain values of the ratio J/U such that the system is su- 
perconducting (or in a Bose glass state) for small Ayly , 
but insulating for large AAij. As discussed earlier, an in- 
creasing value of AAij can be interpreted as an increasing 
value of magnetic field f'^i^/a? for a shghtly disordered 
Josephson junction array in a transverse magnetic field 
equal, on average, to an integer number / of flux quanta 
per plaquette. Thus, our results suggest that, for certain 



values of J/U and integer /, the system undergoes a su- 
perconductor (or Bose glass) to insulator transition as / 
increases. Since a given array would be expected to have 
a fixed value of J /U , such an array may undergo an S-1 
(or BG-1) transition as a function of integer / if J/U is in 
the appropriate range. Our results may not be directly 
applicable to a realistic thin superconducting film in a 
magnetic field because such a film is unlikely to have the 
topology of a Josephson junction network. However, the 
two could exhibit similar phase diagrams. 

In Fig. [THl we have shown a possible phase diagram for 
this system, based on all the numerical data we have ac- 
cumulated. We have drawn the diagram to suggest that 
the entire ordered region for AAij 7^ is of the Bose 
glass type, rather than the S type with z = I. This 
point is discussed further in the next section. Despite 
this assumption, we have, for AAij < 0.854, obtained 
Kc from our calculations of the helicity modulus, as dis- 
cussed above. For AAij — 00, we have used our scaling 
calculations based on the glass order parameter g. For 
AAij = 1-0, we used both methods. They give slightly 
different values of Kc at the phase boundary. It is con- 
ceivable (but, we believe, unlikely) that there is another 
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FIG. 18: Calculated inverse critical coupling constant 1/Kc 
as a function of AAij. The filled points denote the calculated 
points, and the dashed line connecting them is a freehand 
interpolation of the data. We denote the entire ordered re- 
gion for AAij ^ as "BG," consistent with what we believe 
to be the most probable nature of the ordered state. For 
AAij < 0.854, the data come from calculations of the helic- 
ity modulus, as described in the text; for AAij = oo, they 
come from calculation of the glass order parameter g, and for 
AAij = 1.0, they come from both, as shown in the Figure. 



phase boundary separating the S and BG regions some- 
where around AAij = 1-0. Our reasons for beUeving this 
scenario to be unhkely are given in the discussion befow. 



E. Conductivity 

If the transition in our model is from a Mott insulator 
(I) to a superconductor (S), then the helicity modulus 7 
is finite in the S state but vanishes in the state I. Pre- 
cisely at the transition, 7 becomes linear in frequency, 
and the conductivity at the transition can be extracted 
by a scaling analysis [ij , as we review below. In what fol- 
lows, we carry out the scaling analysis over the full range 
of AAij, whether the ordered state is S or BG. 

In order to obtain the value of the conductivity at the 
transition, we need the generalization of the scaling for- 
mulas to frequency-dependent 7. When there is such a 
frequency dependence, Eq. ((28|) is generahzcd to Q 



KM-fik) = GiM^/^S, kM), 



(33) 



where k — 2Trn/M and n is an integer. Precisely at if = 
Kc, G will be a function of only kM , since K — Kc = 0. 
Thus we can introduce another scaling function P, in 
terms of which Eq. ([33|) can be simplified to 



KM-fik) = P{kM). 



From Eq. ((23|) , the conductivity is obtained by taking the 
limit fc ^ after first taking the limit M —^ 00 with a 
small k 1], so that P{kM) ~ kM in the fimit M 00. 
Using the scaling function P, Eq. can be rewritten 



a* PjkM) 
— ~ In fim — m~- 

(JQ kM~^oo kM 



(35) 



where we have also used the relation p{k) = Kj{k). Since 
this quantity is to be calculated for k ^ after M ^ 00, 
the ratio a* jaq will be finite only if the scaling function 
P{x) oc X in this regime. Since k = 2TTn/M , the scaling 
form ([55)1 can be written again as [l| 



cr(n) P(27rn) 



(36) 



This scaling form is expected to be valid only in the 
regime 1 <^ n <^ M . Since it is difficult to carry out cal- 
culations for M large enough that these inequalities are 
satisfied, especially for a disordered system, it is neces- 
sary to incorporate corrections to scaling and express a 
in terms of n and M separately. Since the corrections to 
scaling vanish in the limit n/M — > 00, we expand cr(n) as 
a function of n and M/n using the same form assumed 
in Ref. [l|, namely 



a{n,M/rL) a* 



(37) 



where d and a are fitting constants. The universal con- 
ductivity cr* is found by plotting a{n, M/n) as a function 
of the scaling variable (a/n — n/M) for several lattice 
sizes M and finding the optimal value of a which pro- 
duces the best data collapse onto a single curve. The 
universal conductivity for this value of AAij is the value 
of cr* at which a/n — n/M = 0. 

Using this method, we find the following universal 
conductivities for different values of AAij-. a* /oq — 
0.340 ± 0.006 when AA^ = 1/2, g* /aq = 0.560 ± 0.009 
when A^,j = 1/V2, g* /uq = 1.141±0.088 when AA^ = 
0.854, and c7*/aQ = 1.055 ± 0.090 when AA,j = 00. At 
each of these values of AAij, we apply the method just 
described to calculate the universal conductivity at the 
corresponding Kc values obtained earlier. The results are 
shown in Figs. [Tl HOI EB and[21 respectively. The opti- 
mal values of a's which yield these universal conductiv- 
ities are a = 0.55, 0.19, 0.06, and 0.01 for AAij = 1/2, 
AAij = I/V2, AAij = 0.854, and AAij = 00, respec- 
tively. The accuracy of the calculated o* /gq becomes 
progressively worse as A^l^ increases. In fact, we need 
to obtain the results for AAij — 0.854 and AAij = 00 by 
extrapolation of a(n,M /n) / gq to the optimal values of 
a, using Figs. [2T] and [22l 

The critical coupling constant increases monoton- 
ically with increasing AAij, as we have already noted. 
Similarly, the universal conductivity g* also appears to 
(34) increase monotonically with AAij ■ From these universal 
conductivities, we can plot the universal resistivities as a 



function of AAij ■ These are shown in Fig. [23l The resis- 
tivity decreases with increasing AAij a-U along the phase 
boundary between the Mott insulator and the phase-or- 
dered state. The points represent the results obtained 
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FIG. 19: The conductivity (T(n, M/n) divided by ctq as a 
function of the variable a/n — n/M for several lattice sizes 
M when Ay4ij = 1/2 and /C = /Cc = 0.491. The optimal 
a used here is 0.55. The universal conductivity a* /gq is 
given by that value of a{n, M /n) for which a/n — n/M — 
0, as indicated by the vertical dashed line. The universal 
conductivity thus obtained is cr* /ctq = 0.340 ± 0.006. 




oc/n - n/M 



FIG. 20: Same as Fig. \W\ except that AA.j = 1/^2 and 
Kc = 0.533. The optimal a used here is 0.19. The universal 
conductivity is a'/aq = 0.560 ± 0.009. 



from the QMC simulations, while the dashed lines repre- 
sent a guide to the eye. 

As noted earlier, increasing AAij in our model corre- 
sponds to increasing magnetic field for a shghtly disor- 
dered Josephson junction array in a uniform transverse 
magnetic field. Thus, our results which show a decrease 
in the universal resistivity with increasing AA^, prob- 
ably cannot be directly compared to the experiments 
reported in Refs. [1^ [2lj and the numerical results of 
Ref. p^ . The experiments consider a disordered Bi film 
in a uniform magnetic field rather than a slightly disor- 
dered Josephson junction array in a periodic field. How- 
ever, both the experiments and our calculations find a 
zero-temperature magnetic-field-tuned transition from a 
phase-ordered state to an insulator, and both the exper- 
imental papers and previous calculations interpret this 
phase transition using scaling fits of the low-temperature 
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FIG. 21: Same as Fig. [191 except that AAij = 0.854 and Kc = 
0.585. The optimal a used here is 0.06. We therefore have 
to extrapolate the plot of the conductivity a{n, M /n) / gq to 
reach the point at a/n— n/M = 0. The universal conductivity 
thus obtained is a* /aq = 1.141 ± 0.088. 




FIG. 22: Same as Fig. EH except that AAij =oo,Kc = 0.630, 
and z — 1.1. The optimal a used here is 0.01. The universal 
conductivity is also obtained from extrapolation of the data, 
resulting in a'/ag = 1.055 ± 0.090. 

transport properties near the critical field, as we do here 
for our model. 



IV. DISCUSSION 

In this paper, we have calculated the transition be- 
tween the superconducting and the insulating state for a 
model disordered 2D superconductor in a magnetic field. 
We treat the superconductor as a square Josephson junc- 
tion array with an intergranular Josephson coupling en- 
ergy J and a finite capacitive energy described by an 
on-site charging energy U. To include the effect of both 
disorder and a transverse magnetic field, we include a 
random magnetic phase factor Aij in the Josephson cou- 
pling between grains; Aij is assumed to have a Gaus- 
sian distribution with zero mean and a root-mean-square 
width AAi-j. 
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FIG. 23: The universal resistivity p* divided by pq as a func- 
tion of AAij. In each case p* /pQ = o-g/cr* . The dashed lines 
are cubic spline fits to the data. "I" and "BG" denote the 
insulating and Bose glass phases, respectively. The supercon- 
ducting phase, in our interpretation, occurs only at AAij — 0. 
The open circle at AAij = 1.0 is an average between the two 
values of p*/pQ from the two dashed lines as in Fig. [T51 



Although our model is certainly artificial, it should rea- 
sonably represent the effects of a magnetic field applied 
to a disordered 2D array at integer /. Specifically, the 
model resembles a spatially disordered 2D granular su- 
perconductor in a uniform magnetic field, in which the 
plaquettes have slightly different areas. For small AAij, 
the root-mean-square frustration per plaquette is small, 
but at all nonzero AAij , the grain plaquettes are ran- 
domly frustrated, a feature which should be relevant in 
real 2D films. Our model provides a way of interpolat- 
ing smoothly between the zero-field and high-field limits 
[62 |. We have confirmed numerically that, at least in the 
high-field limit, the two models give a BG-I transition at 
the same value of K. 

Our numerical results suggest that, for any value of 
AAij , the system undergoes a transition from an insu- 
lating (I) state to an ordered state. We believe that the 
transition is I to Bose glass (BG) over the entire range of 
AAij, except for AAij = 0. Supporting this hypothesis 
is the fact that g exhibits excellent scaling at AAij — oo 
and very good scaling at AAij — 1-0. This hypothesis 
is also indirectly supported by the fact that the scaling 
behavior of 7 becomes progressively worse as AAij in- 
creases. To provide stronger support for this hypothesis 
numerically for smaller AAij, we would need to go to 
much larger Monte Carlo sample sizes, using the correct 
value of 2; ^ 1.3. 

In support of this scenario, we now describe a sim- 
ple argument, based on the well-known Harris criterion 
[63] , which suggests that the S state is unstable against a 
small random perturbation of the type we consider here. 
Consider the zero-field version of model ([T]), in the pres- 
ence of some kind of weak uncorrelated disorder. In a 
region of size ^'^ (where d = 2 and ^ is the correlation 
length of the unperturbed system), the critical value of 



K should fiuctuate by an amount of order . Near 
the critical value of K, the correlation length of the un- 
perturbed system varies with K according to the relation 
f oc \{K — Kc) / K^~'^ . In order for the transition of the 
unperturbed system to be unaffected by the disorder, the 
Harris criterion suggests that v > 2/d. For the present 
case, d = 2, but because we are dealing with a quan- 
tum transition, v is that of a 3D XY transition, namely 
1/^2/3 [131 ■ Thus, the inequality is not satisfied and 
we expect this quantum phase transition to be unstable 
against point disorder in 2D. 

For the present model, the randomness is indeed un- 
correlated within the plane, as required by the above ar- 
gument, but it is somewhat different from the usual point 
disorder. For small AAij, the random part of the Hamil- 
tonian may be written 

K^^iAA,) ~ - sin(0, - 0j)AA,, (38) 

where AAij is a Gaussian random variable. Despite the 
form of this disorder, it seems reasonable that the disor- 
der would have at least as strong an effect on the phase 
transition of the pure model as more conventional point 
disorder. Therefore, we suggest, based on this rough ar- 
gument, that the 3D XY phase transition of the pure 
model is unstable against this random field perturbation 
for arbitrarily weak AAij ■ 

The next question is, to what is the 3D XY tran- 
sition unstable? The most likely scenario is that the 
transition is of the I to BG class over the entire range 
AAij 7^ 0. The seeming presence of the S phase at 
small but finite AAij is probably due to the fact that 
our samples are not large enough to exhibit the expected 
BG phase. For example, at AAij = 1/2, the rms vari- 
ation in total flux through a single plaquette would be 
2 • (l/2)$o/(27r) = $o/(27r). Thus, the rms variation 
in total flux through a lattice of plaquettes would be 
L$o/(27r) (where the factor of L comes from the fact that 
4L is the perimeter of the plaquettes), and hence, even 
for L = 12, would be only about two flux quanta. This 
value is not large enough to yield results characteristic of 
AAij = 00, at which we have shown the Bose glass is the 
stable ordered phase. Thus, indeed, the sample sizes we 
have considered are simply not large enough to exhibit 
fully developed Bose-glass scaling at AAij ~ 1/2, or even 
at larger values than this. Nonetheless, we have the basic 
result that, for any AAij, there is a transition from an in- 
sulating state to an ordered state with increasing values 
of the coupling parameter K = ^ J /{2U). As the above 
argument suggests, we believe that this ordered state is 
a BG phase for any nonzero AAij. 

The critical coupling constant for the transition 
from I to BG increases monotonically with increasing 
AAij. Thus, for certain values of K, the material is 
in the BG state at low AAij, but goes through a BG- 
I transition as AAij increases. In a disordered material, 
increasing AAij can be identified with increasing trans- 
verse magnetic field for a slightly disordered Josephson 
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junction array at integer /, in the sense which we have 
discussed earher. Since any given material should have 
a fixed if, independent of AAij , this trend implies that 
some materials, which have suitable values of K, will go 
through a BG-I transition with increasing field. A mate- 
rial with a smaller K will remain insulating for all AAij , 
while one with a larger K, would remain a BG for all 
fields. All this behavior follows from the phase diagram 
drawn in Fig. [181 In each case, the I phase in our model 
is a Mott insulator, since Cooper pairs are localized by 
Coulomb repulsion rather than by disorder. 

Besides calculating the critical values of coupling con- 
stant Kc, we have also computed the universal conductiv- 
ities a* as a function of AAij for these transitions. The 
values of both Kc and p* / PQ = (o'*/ctq)~^ are shown 
for various values of AAij in Fig- [131 In all cases, these 
values are obtained by a scaling analysis of the numeri- 
cally calculated helicity modulus 7 and the renormalized 
coupling constant g. 

Our results maybe consistent with experimental find- 
ings as in Refs. P, HU, [H, HI jH, [H, U 113, H . 
The data in these references indicate that a-InO^ films 
[H, [111, granular In films a-MoGe films Bi 

films~p,l2l|, Nd2-xCe^Cu04-tj, films [IH, TiN films 
[3^. l34| ■ and Nbo.15Sio.85 films [sl] show that the resis- 
tance per square, normalized by the resistance at the 
transition, at very low temperatures decreases as a func- 
tion of the scaled magnetic field B when the magnetic 
field is less than the critical value Be, while it increases 
when B > Be- It should be kept in mind, of course, that 
our model calculations refer to a disordered Josephson 
array at an integer number of flux quanta per plaquette, 
on average, whereas the experiments deal with systems 
having a possibly different topology. 

Numerically, there are several ways in which our cal- 
culations could be further improved. In some cases, the 
number of realizations (100) we have used for Aij may 
be insufficient to provide accurate statistics and might 
lead to significant numerical uncertainties. Our choice 
for the number of realizations is dictated by a compro- 
mise between computing costs and statistical errors. Be- 
cause of the large amount of computing time involved, 
we have carried out our calculations only up to a lattice 



size at most of 20 on an edge, and have considered only 
five different nonzero AA^'s. Our results would have had 
greater accuracy and given a more detailed picture of the 
phase diagram if we had been able to include more values 
of AAij, a larger number of realizations, and, especially, 
larger lattice sizes. In addition, our "world-line" algo- 
rithm [33,IH,|6^ could be replaced by other approaches, 
such as a "worm" algorithm [s^l or a stochastic series 
expansion [III 113, HH , possibly leading to better con- 
vergence. It might also be valuable to develop another 
model in which the disorder is introduced in a manner 
closely resembling that in actual superconducting ffims. 
Finally, we note that the same approach could be used 
to calculate the finite-frequency conductivity of the low- 
temperature phase for various values oiAA.j ^M,^- 
To summarize, we have carried out extensive quantum 
Monte Carlo simulations of a model for a transition from 
a Mott insulator to a superconducting phase at low tem- 
peratures. The model is characterized by a continuously 
tunable disorder parameter AAij ■ Our numerical results 
suggest that, for any nonzero AAij, there is such a tran- 
sition, and that the ordered phase is a Bose glass. The 
evidence that the ordered phase is a Bose glass is strong 
for AAij = 00, but less conclusive for smaller AAij. We 
also find that, for certain values of the coupling variable 
K, the system can go from BG at small AAij to a Mott 
insulator at large AAij. We also discuss the possibility 
that this transition may be related to the field-driven su- 
perconductor to insulator transition seen in a number of 
superconducting films. It would be of great interest if 
our results could be compared to a suitable experimental 
realization. 
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